function [out,excessHoldPerReturn,termSpreads] = forecastRegressionYieldSpread(yields,shortRateForXhr,shortRateOneYear,maturities,REC,h)
% Data inputs:
% yields: T * numYields where yields are annually and continuously compounded
% shortRateForXhr: Short rate for excess returns
% maturities: 1*numYields with maturities in months
% RECdummy: 0 1 regression dummy, where 0 is boom and 1 is recression
% h: sampletime_per_maturityTime

% Preable
T         = size(yields,1);
numYields = length(maturities);

% Term spread: y[t](n)-y[t](1-year) for all n.
termSpreads      = yields(:,1:end)-repmat(shortRateOneYear,1,numYields);

% Ex post holding period returns: hr[t]_n = n*y[t-h](n) - (n-h)*y[t](n-h)
annualMaturities = maturities/12;
holdPerReturn   =[nan(T-h,1), (repmat(annualMaturities(1,2:end),T-h,1).*yields(1:end-h,2:end)...
    -repmat(annualMaturities(1,1:end-1),T-h,1).*yields(1+h:end,1:end-1))];

% Ex post excess holding period returns: xhr[t]_n = (hr[t]_t - y[t-h](h))*1/h           
excessHoldPerReturn = (holdPerReturn/h - repmat(shortRateForXhr(1:end-h,1)/12,1,numYields));               
               
%% Runing the forecast regressions for all maturities and for the three predictors
%index_RECh       = find(REC(1:end-h,1)==1);
%index_EXPh       = find((1-REC(1:end-h,1))==1);
beta             = NaN(2,numYields-1);
beta_tstat       = NaN(2,numYields-1);
R2               = NaN(1,numYields-1);
betaRegime       = NaN(4,numYields-1);
betaRegime_tstat = NaN(4,numYields-1);
R2regime         = NaN(1,numYields-1);
startMat         = find(maturities >= 2*12);

%for n=1:numYields-1
for n=startMat:numYields
    % excessHoldPerReturn regressed on term spread
    Y = excessHoldPerReturn(:,n);
    
    X = [ones(size(Y,1),1) termSpreads(1:end-h,n)];
    
    res            = nwest(Y,X,2*h);
    beta(:,n)      = res.beta;
    beta_tstat(:,n)= res.tstat;
    R2(1,n)        = res.rsqr;
    %residuals      = Y - X*res.beta;
    %residualsConst = Y - mean(Y);
    
    % Regime dependent: excessHoldPerReturn regressed on term spread
    %X1=[X.*repmat(1-REC(1:end-h,1),1,2) X.*repmat(REC(1:end-h,1),1,2) ];
    X1=[ones(size(Y,1),1) REC(1:end-h,1) termSpreads(1:end-h,n) termSpreads(1:end-h,n).*REC(1:end-h,1) ];
    resRegime             = nwest(Y,X1,2*h);
    betaRegime(:,n)       = resRegime.beta;
    betaRegime_tstat(:,n) = resRegime.tstat;
    R2regime(1,n)         = resRegime.rsqr;
    
    %R2_rec(1,n) = (1-sum(residuals(index_RECh).^2)/sum(residualsConst(index_RECh).^2))*100;
    %R2_exp(1,n) = (1-sum(residuals(index_EXPh).^2)/sum(residualsConst(index_EXPh).^2))*100;
end

%% Output
out.beta             = beta;
out.beta_tstat       = beta_tstat;
out.beta_se          = beta./beta_tstat;
out.R2               = R2;
out.betaRegime       = betaRegime;
out.betaRegime_tstat = betaRegime_tstat;
out.betaRegime_se    = betaRegime./betaRegime_tstat;
out.R2regime         = R2regime;
out.maturities       = maturities;


end

